#ifndef UHM_DISCRETE_BASE_MESH_HXX_
#define UHM_DISCRETE_BASE_MESH_HXX_

#define UHM_MESH_ELEMENT_RESERVE_SIZE 100

namespace uhm {

  class Mesh_Base_ {
  public:

    /*
    enum { RELATION_ME2PAR=1,
           RELATION_PAR2ME } Relation_Direction_;
    */

    // Be careful.. uhmfile think default is 0
    enum { NODE_DEFAULT=0,
           NODE_BOUNDARY } Kind_;

    enum { NODE_NOT_SEPARATED=0,
           NODE_FACTOR,
           NODE_SCHUR } Node_Status_;

    enum { ELEMENT_MARKER_ZERO=0,
           ELEMENT_REUSED } Element_Marker_;

    Mesh_Base_() { }
    Mesh_Base_(const Mesh_Base_ &b) { }
    virtual~Mesh_Base_() { }
  };
  
  template<class Id_Nd_, class Id_El_, class Uhm_>
  class Mesh_ : public Mesh_Base_,
                public Disp_ {

  public:    

    typedef Id_Nd_ _Id_Nd_;
    typedef Id_El_ _Id_El_;
    typedef Uhm_ _Uhm_;

    // Forward declaration
    class Relation_;
    class Node_;
    class Element_;
    
    class Relation_ : public Disp_ {
    public:
      
    protected:
      Mapper_<Id_Nd_> _me2par_s2f, _me2par_s2s;
      Int_ flag;
    public:
      Relation_() :  _me2par_s2f(), 
                     _me2par_s2s() { flag = 0; }
      virtual~Relation_() {  }

      void set(Element_ *me, Element_ *par) {
        _me2par_s2f.set(&me->schur(),  &par->factor());
        _me2par_s2f.relate();

        _me2par_s2s.set(&me->schur(),  &par->schur());
        _me2par_s2s.relate();

        flag = 1;
      }

      void clear() {
        _me2par_s2f.clear();
        _me2par_s2s.clear();
      }
      
      void get(Mapper_<Id_Nd_> &f, Mapper_<Id_Nd_> &s) {
        f = _me2par_s2f; 
        s = _me2par_s2s; 
      }

      virtual std::ostream& content(std::ostream& os) const  {
        if (!flag) return os;
        return (os 
                << "- Relation_ -\n" 
                << " me2par (s2f) = \n" << _me2par_s2f << "\n"
                << " me2par (s2s) = \n" << _me2par_s2s << "\n");
      }
    };
    
    // Node
    // ----
    class Node_ : public Disp_ {
    private:
      // Id include dof and misc info such as p.
      Id_Nd_ _id;
      Int_ _kind;

      std::set< Element_* > _owner;

    protected:
      Int_ _offset, _marker; 

      // offset - global assembly offset
      // marker - temporary marker 

    public:

      Node_() :  _kind(NODE_DEFAULT), 
                 _offset(-1), 
                 _marker(0) { }
      Node_(Id_Nd_ id) : _kind(NODE_DEFAULT), 
                         _offset(-1), 
                         _marker(0), 
                         _id(id) { }
      Node_(Id_Nd_ id, Int_ kind) : _kind(kind), 
                                    _offset(-1),
                                    _marker(0), 
                                    _id(id) { }
    
      Node_(const Node_ &b) : _kind(b._kind),
                              _offset(b._offset),
                              _marker(b._marker),
                              _id(b._id),
                              _owner(b._owner) { }
      virtual~Node_() { }

      Int_ set_offset(Int_ offset) { 
        _offset = offset; 
        return (offset + _id.dof());
      }
      Int_ get_offset() const { return _offset; }
      Int_ offset() { return _offset; }
      
      void set_kind(Int_ kind) { _kind = kind; }
      Int_ get_kind() const { return _kind; }
      Int_ kind() { return _kind; }
      
      void set_marker(Int_ marker) { _marker = marker; }
      Int_ get_marker() const { return _marker; }
      Int_ marker() { return _marker; }
      
      Id_Nd_ get_id() const { return _id; }
      Id_Nd_& id() { return _id; }

      std::set< Element_* >& owner() { return _owner; }
      
      Int_ get_n_owner() const { return _owner.size(); }
      void add_owner(Element_* e) {
#pragma omp critical 
        {
          _owner.insert(e);
        }
      }
      void remove_owner(Element_* e) {
#pragma omp critical
        {
          _owner.erase( _owner.find(e) );
        }
      }

      void clear() {
        _marker = 0; _offset = 0; _owner.clear();
      }
      virtual Bool_ operator<(const Node_ &b) const  {
        //Int_ val = this->_kind - b._kind;
        //return (val ? val < 0 : this->_id < b._id);
        return (this->_id < b._id);
      }
      virtual Bool_ operator==(const Node_ &b) const {
        //return (this->_id == b._id && this->_kind == b._kind );
        return (this->_id == b._id);
      }
      virtual Bool_ operator!=(const Node_ &b) const  {
        return !(*this == b);
      }
      virtual std::ostream& content(std::ostream& os) const  {
        os << "- Node_ " 
           <<  _id 
           << ", Kind(" << _kind << ")"
           << ", Marker(" << _marker << ")"
           << ", Offset(" << _offset << ")";

        if (_owner.size()) {
          os << ", Owners ("<< _owner.size() << ") ";
          typename std::set< Element_* >::iterator it;
          for (it=_owner.begin();it!=_owner.end();++it)
            os << "[ " << (*it)->get_id() << " ]";
        }
        return (os);
      }
    };
    
    class Element_ : public Disp_ {
    public: 
      class Entry_ : public Disp_ {
        Node_* _node;
        Int_   _status;
      public:
        Entry_() : _node(NULL), _status(NODE_NOT_SEPARATED) { }
        Entry_(Node_* n, Int_ status) : _node(n), _status(status) { }
        Entry_(const Entry_ &b) : _node(b._node), _status(b._status) { }

        void set_node(Node_* n) { _node = n; }
        Node_* get_node() const { return _node; }
        Node_* node() { return _node; }

        void set_status(Int_ status) { _status = status; }
        Int_ set_status() {
          UHM_ERROR(_node->get_n_owner(), UHM_ERROR_MESH_NODE_SET_STATUS);
          if (_node->get_n_owner() > 1 ||
              _node->get_kind() == NODE_BOUNDARY) 
            _status = NODE_SCHUR;
          else
            _status = NODE_FACTOR;
          return _status;
        }
        Int_ get_status() const { return _status; }
        Int_ status() { return _status; }

        Entry_& operator=(const Entry_ &b) {
          _node = b._node;
          _status = b._status;
          return *this;
        }
        Bool_ operator<(const Entry_ &b) const {
          return (_node->get_id() < b._node->get_id());
        }
        Bool_ operator==(const Entry_ &b) const {
          return (this->_node->get_id() == b._node->get_id());
        }
        Bool_ operator!=(const Entry_ &b) const  {
          return !(*this == b);
        }

        virtual std::ostream& content(std::ostream& os) const  {
          return (os << "Status:(" << _status << ") "<< *_node) ;
          //return (os << "(" << _node->get_id() << ", Status(" << _status << ")");
        }
      };

    private:
      Id_El_ _id;
      std::vector< Entry_ > _nds;
      Map_< Id_Nd_ > _factor, _schur;

      // For leaf, relation explains bijection to work storage.
      // Non leaf will related to its upper and lower hierarchy.
      Relation_ _relation;

      Uhm_ _uhm;

      // recorded flop for statistics. if uhm is reused, _flop = 0.0
      Double_ _flop;

    protected:
      Int_ _marker, _group;

    public:
      Element_() : _marker(0), 
                   _group(0), 
                   _flop(0.0) { }
      Element_(Id_El_ id) : _id(id), 
                            _marker(0), 
                            _group(0),
                            _flop(0.0) { }

      // Watch out that I don't copy relation since it is composed of 
      // pointer objects. It is ambiguous if they are copied all over...
      Element_(const Element_ &b) : _id(b._id), 
                                    _nds(b._nds), 
                                    _factor(b._factor), 
                                    _schur(b._schur), 
                                    _uhm(b._uhm),
                                    _marker(b._marker),
                                    _group(b._group),
                                    _flop(b._flop) { }

      void set_marker(Int_ marker) { _marker = marker; }
      Int_ get_marker() const { return _marker; }
      Int_& marker() { return _marker; }

      void set_group(Int_ group) { _group = group; }
      Int_ get_group() const { return _group; }
      Int_& group() { return _group; }

      void set_flop(Double_ flop) { _flop = flop; }
      Double_ get_flop() const { return _flop; }
      Double_& flop() { return _flop; }

      Id_El_ get_id() const { return _id; }
      Id_El_& id() { return _id; }

      Uhm_ get_uhm() const { return _uhm; }
      Uhm_& uhm() { return _uhm; }

      // Manual access
      Map_<Id_Nd_>& factor() { return _factor; }
      Map_<Id_Nd_>& schur() { return _schur; }
      std::vector<Entry_>& nds() { return _nds; }

      Relation_& relation() { return _relation; }

      void nds2map(Map_<Id_Nd_> &factor, Map_<Id_Nd_> &schur) {
        factor.push_begin(0);
        factor.push_end();

        schur.push_begin(_nds.size());
        for (Int_ i=0;i<_nds.size();++i) 
          schur.push(_nds[i].node()->id());
        schur.push_end();
      }



      Bool_ is_separated() const { 
        return (_nds.size() == (_factor.get_size() + _schur.get_size())); 
      }

      void separate() {
        // reset marker
        this->set_marker(ELEMENT_MARKER_ZERO);

        Int_ cnt = 0;
        typename std::vector<Entry_>::iterator it;
        for (it=_nds.begin();it!=_nds.end();++it) {

          // set status
          Node_* n = it->node();
          if (n->get_n_owner() == 1 &&
              n->get_kind() != NODE_BOUNDARY ) 
            it->set_status(NODE_FACTOR);
          else 
            it->set_status(NODE_SCHUR);
          
          // if it is factor add cnt
          cnt += (it->get_status() == NODE_FACTOR);
        }

        _factor.clear(); _factor.set_size(cnt);
        _schur.clear();  _schur.set_size(_nds.size()-cnt);
        
        _factor.push_begin();
        _schur.push_begin();
        for (it=_nds.begin();it!=_nds.end();++it) {
          Node_* n = it->node();
          switch (it->get_status()) {
          case NODE_FACTOR: _factor.push(n->id()); break;
          case NODE_SCHUR:  _schur.push(n->id());  break;
          case NODE_NOT_SEPARATED:
            UHM_ERROR(false, UHM_ERROR_MESH_ELEMENT_SEPARATE);
          }
        }
        _factor.push_end();
        _schur.push_end();
      }

      void leave_boundary_nodes_only() {
      }

      // This assume that element look down to children
      // Not parent, if this element update parent, it will reset every time
      // child update the parent.
      void push_begin(Int_ size=UHM_MESH_ELEMENT_RESERVE_SIZE) {
        _nds.clear(); _nds.reserve(size);
      }
      void push(Node_* n) {
			_nds.push_back(Entry_(n, NODE_NOT_SEPARATED));
      }
      void push_end() {
        std::sort(_nds.begin(), _nds.end());
        typename std::vector< Entry_ >::iterator it;
        it = std::unique(_nds.begin(), _nds.end());
        _nds.resize( it - _nds.begin() );
        
        for (it=_nds.begin();it!=_nds.end();++it)
          it->node()->add_owner(this);
      }

      // update element from all children, should be wrapped with push_begin, push_end.
      void update(Element_ *c) {
        UHM_ERROR(c->is_separated(), UHM_ERROR_MESH_ELEMENT_UPDATE);

        std::vector<Entry_> &c_nds = c->nds();
        for (Int_ i=0;i<c_nds.size();++i) {
          Node_* n = c_nds[i].node();
          if (c_nds[i].get_status() == NODE_SCHUR) {
            n->remove_owner(c);
            this->push(n);
          }
        }
      }

      Int_ set_nodes_offset(Int_ offset) {
        // Set offsets for factoring nodes.
        typename std::vector< Entry_ >::iterator it;
        for (it=_nds.begin();it!=_nds.end();++it)
          offset = it->node()->set_offset(offset);
        return offset;
      }

      void clear(Int_ is_leaf) {
        _marker = 0; _flop =0.0;
        _relation.clear();
        _factor.clear();
        _schur.clear();
        //_uhm.clear();
        
        if (is_leaf) {
          // if leaf, recover origianl connectivity
          typename std::vector< Entry_ >::iterator it;
          for (it=_nds.begin();it!=_nds.end();++it) {
            it->set_status(NODE_NOT_SEPARATED);
            it->node()->add_owner(this);
          }
        } else {
          // otherwise, clear all
          _nds.clear();
        }
      }

      virtual Bool_ operator<(const Element_ &b) const  { return (this->_id < b._id); }
      virtual Bool_ operator==(const Element_ &b) const { return (this->_id == b._id);}
      virtual Bool_ operator!=(const Element_ &b) const { return !(*this == b); }

      virtual std::ostream& content(std::ostream& os) const  {
        os << "- Element_ " << _id 
           << ", Marker(" << _marker << ")"
           << ", Group(" << _group << ")";
        
        if (_nds.size()) {
          os << "\n- Nodes ("<< _nds.size() << ") -\n";
          for (Int_ i=0;i<_nds.size();++i)
            os << _nds.at(i) << std::endl;
        }
        if (_factor.get_size())
          os << "\nFactor " << _factor;
        if (_schur.get_size())
          os << "\nSchur  " << _schur ;

        os << _relation;
        // os << "\nMatrices " << _uhm;

        return os;
      }

    };
  private:
    // This id goes negative in order not to have duplication with FEM interface.
    Id_El_ _el_id;

  protected:
    std::set< Node_ > _nodes;
    std::set< Element_ > _elements;
    
  public:

    Mesh_() : Mesh_Base_() { }
    Mesh_(Id_El_ el_id) : Mesh_Base_(), _el_id(el_id) { }
    Mesh_(const Mesh_ &b) : Mesh_Base_(b), _el_id(b._el_id) { }
    
    void clear() {
      _el_id = Id_El_();
      _nodes.clear();
      _elements.clear();
    }

    Node_* add_node(Id_Nd_ id, Int_ kind=NODE_DEFAULT) {
      std::pair< typename std::set< Node_ >::iterator, Bool_ > ret;      

#pragma omp critical
      {
        ret = _nodes.insert( Node_(id,kind) ) ;
      }

      return &((Node_&)*ret.first);
    }

    // Id does not compare dof
    Node_* find_node(Id_Nd_ id) {
      typename std::set< Node_ >::iterator it;
      it = _nodes.find(Node_(id)) ;

      return (it == _nodes.end() ? NULL : &((Node_&)*it));
    }

    void clear_nodes() {
      typename std::set< Node_ >::iterator it;
      for (it=_nodes.begin();it!=_nodes.end();++it) {
        Node_* n = (Node_*)&(*it);
        n->clear();
      }
    }
    
    Element_* add_element() {
      std::pair< typename std::set< Element_ >::iterator,Bool_ > ret;
#pragma omp critical
      {
        --_el_id;
        ret = _elements.insert( Element_(_el_id) );
        UHM_ERROR( ret.second, UHM_ERROR_MESH_ADD_ELEMENT);
      }
      return &((Element_&)*ret.first);
    }

    Element_* add_element(Id_El_ id) {

      std::pair< typename std::set< Element_ >::iterator, Bool_ > ret;
#pragma omp critical
      {
        ret = _elements.insert( Element_(id) );
        UHM_ERROR( ret.second, UHM_ERROR_MESH_ADD_ELEMENT);
      }
      return &((Element_&)*ret.first);
    }

    Element_* find_element(Id_El_ id) {
      typename std::set< Element_ >::iterator it;
      it = _elements.find(Element_(id)) ;

      return (it == _elements.end() ? NULL : &((Element_&)*it));
    }


    Int_ get_n_nodes() const { return _nodes.size(); }
    Int_ get_n_elements() const { return _elements.size(); }
    Int_ get_n_elements_created_by_user(Id_El_ mesh_define_begin) const {
      Int_ cnt=0;
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        if (e->id() > mesh_define_begin)
          ++cnt;
      }
      return cnt;
    }

    void clear(Id_El_ el_id_begin) {
      _el_id = el_id_begin;
      _nodes.clear();
      _elements.clear();
    }

    void dump(std::vector<Element_*> &els, Int_ pid=0) {
      els.clear();
      els.reserve(this->_elements.size());
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        if (e->id().get_part() == pid)
          els.push_back(e);
      }
    }
    
    void dump(std::vector<Node_*> &nds, Int_ pid=0) {
      nds.clear();
      nds.reserve(this->_nodes.size());

      typename std::set< Node_ >::iterator it;
      for (it=_nodes.begin();it!=_nodes.end();++it) {
        Node_* n = (Node_*)&(*it);
        nds.push_back(n);
      }
    }

    void set_markers(Int_ marker) {
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        e->set_marker(marker);
      }
    }

    // Not involve with pid
    // need to think what needs pid
    void set_matrices(Int_ type, Int_ rhs, Int_ blk) {
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        e->uhm().set(type, e->factor().dof(), e->schur().dof(), rhs, blk);
      }
    }

    // This does not create ABR although element is a leaf
    // Should be managed in merging
    void create_matrices() {
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        e->uhm().create_flat_without_buffer();
        e->uhm().create_flat_buffer(); // do not chreate ABR
        e->uhm().create_hier();
      }
    }
    
    void free_matrices() {
      typename std::set< Element_ >::iterator it;
      for (it=_elements.begin();it!=_elements.end();++it) {
        Element_* e = (Element_*)&(*it);
        e->uhm().free_all();
      }
    }

    Int_ get_dof(Int_ &max_dof, Int_ &min_dof) const  {
      Int_ r_val = 0;

      { // ** sum nodal dofs
        Int_ dof = 0;
        typename std::set< Node_ >::iterator it;
        for (it=_nodes.begin();it!=_nodes.end();++it) {
          Node_* n = (Node_*)&(*it);
          dof += n->id().get_dof();
        }
        r_val = dof;
      }

      { // find out max min uhm
        max_dof = 0; min_dof = 10000000;
        typename std::set< Element_ >::iterator it;
        for (it=_elements.begin();it!=_elements.end();++it) {
          Element_* e = (Element_*)&(*it);
          Int_ dof = (e->factor().dof() + e->schur().dof());
          max_dof = (max_dof < dof ? dof : max_dof);
          min_dof = (min_dof > dof ? dof : min_dof);
        }
      }

      return r_val;
    }

    void nodes2map(Map_<Id_Nd_> &nodes) {
      nodes.clear();
      nodes.push_begin(_nodes.size());
      typename std::set< Node_ >::iterator it;
      for (it=_nodes.begin();it!=_nodes.end();++it) {
        Node_* n = (Node_*)&(*it);
        nodes.push(n->id());
      }
      nodes.push_end();
    }

    virtual std::ostream& content(std::ostream& os) const  {
      this->Disp_::content(os);
      os << "- Mesh_ -" << std::endl
         << "  # Nodes    " << get_n_nodes() << std::endl
         << "  # Elements " << get_n_elements() << std::endl;

      {
        typename std::set< Node_ >::iterator it;
        for (it=_nodes.begin();it!=_nodes.end();++it)
          os << (*it) << std::endl;
      }
      os << std::endl;
      {
        typename std::set< Element_ >::iterator it;
        for (it=_elements.begin();it!=_elements.end();++it)
          os << (*it) << std::endl;
      }
      return os;
    }

  };

}


#endif
